## Authors: Alexander Herzog and Kenneth Benoit
## Date: May 30, 2015
## Replication file for JOP article "The Most Unkindest Cuts: Speaker Selection and Expressed Government Dissent During Economic Crisis"

rm(list = ls(all = TRUE))

library(rjags)
set.seed(42)

### PATHS ########################################
dataAll <- "./generated_data/working_data_all.RData"
##################################################

# MCMC values
mcmc.adapt <- 500
mcmc.burn <- 1000
mcmc.sample <- 2000
mcmc.thin <- 2
mcmc.chains <- 2

# JAGS model 
JagsModel <- "model{
	for (i in 1:N)
	{
		y[i] ~ dbern(p[i])
		probit(p[i]) <- a[m[i]] + inprod(b, X[i,1:K])
	}
	
	for (j in 1:J)
	{
		a[j] ~ dnorm(mu.a, tau.a)
	}

	for (k in 1:K)
	{
		b[k] ~ dnorm(0, 0.0001)
	}

	mu.a ~ dnorm(0, 0.0001)
	tau.a <- pow(sigma.a, -2)
	sigma.a ~ dunif(0, 1000)

}"

# function for JAGS data for multilevel models
forJagsData <- function(DV, xMatrix, groupID) {
    list(
        y = DV,
        X = xMatrix,
        N = nrow(xMatrix), 
        K = ncol(xMatrix),
        J = length(unique(groupID)),
        m = groupID
    )
}

load(dataAll)
data <- dataAll


# Model 1 without interaction effects
# -----------------------------------
print("Estimating selection model 1")

eq <- spoke ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis

X <- model.matrix(eq, data)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(data$spoke, X, data$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

spoke1.posterior <- coda.samples(model,
                                 variable.names=c("a", "b", "mu.a", "sigma.a"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(spoke1.posterior,file="./estimated_models/selection_model1_posterior.RData")


# Model 2 with interaction effects
# ---------------------------------
print("Estimating selection model 2")

eq <- spoke ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis + crisis*unemployment + crisis*safetyScaled

X <- model.matrix(eq, data)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(data$spoke, X, data$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

spoke2.posterior <- coda.samples(model,
                                 variable.names=c("a", "b", "mu.a", "sigma.a"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(spoke2.posterior,file="./estimated_models/selection_model2_posterior.RData")


# Model 3 with three-way interaction effects
# ------------------------------------------
print("Estimating selection model 3")

eq <- spoke ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis + crisis*unemployment + crisis*safetyScaled + crisis*government + government*unemployment + government*safetyScaled + crisis*government*unemployment + crisis*government*safetyScaled

X <- model.matrix(eq, data)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(data$spoke, X, data$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

spoke3.posterior <- coda.samples(model,
                                 variable.names=c("a", "b", "mu.a", "sigma.a"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(spoke3.posterior,file="./estimated_models/selection_model3_posterior.RData")
